CRISPR-EATING is a molecular biology protocol to generate libraries of CRISPR guide RNAs. The use of this this approach to generate a library suitable for chromosomal locus imaging requires ways to avoid regions that will be processed into non-specific guides, which (in part) is what these scripts are designed to achieve.
These scripts contain a set of functions that are assembled into a workflow to:
Following the generation of suitable PCR primers from this tool, the "wet" portion of the protocol is as follows:
makeblastdb -in LAEVIS_7.1.repeatmasked.fa -dbtype nucl -parse_seqids -out xl71 -title ‘xl71’
Hsu PD, Scott DA, Weinstein JA, Ran FA, Konermann S, Agarwala V, et al. DNA targeting specificity of RNA-guided Cas9 nucleases. Nat Biotechnol. Nature Publishing Group; 2013;31: 827–832. doi:10.1038/nbt.2647
In [1]:
import Bio
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio import SeqIO
from Bio.Blast import NCBIXML
from Bio import Restriction
from Bio.Restriction import *
from Bio.Alphabet.IUPAC import IUPACAmbiguousDNA
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
import cPickle as pickle
import subprocess
import matplotlib
from eating import *
import multiprocessing as mp
from operator import itemgetter, attrgetter, methodcaller
import numpy
%pylab inline
In [2]:
fasta_file = SeqIO.parse("../../Genomic Data/LAEVIS_7.1.repeatMasked.fa", "fasta")
In [3]:
handle = open("../../Genomic Data/LAEVIS_7.1.repeatMasked.fa", 'rb')
xl71genome = SeqIO.parse(handle, "fasta", alphabet=IUPACAmbiguousDNA())
xl71genomedict = {}
for item in xl71genome:
xl71genomedict[item.id] = item
In [4]:
len(xl71genomedict)
Out[4]:
The FASTA file we've loaded (fasta_file) contains the entire X. laevis genome. The X. laevis genome hasn't yet been definitively assembled into physical chromosomes - instead, it's a large number of contigs or "scaffolds". For the purposes of making a library that labels a single region, we want to work with a big piece that we know is contiguous. So, we find the longest "scaffold".
In [5]:
longest = 0
for item in fasta_file:
if len(item) > longest:
longest = len(item)
longscaffold = [item]
In [6]:
print(longscaffold[0].name + " is the longest scaffold at " "{:,}".format(len(longscaffold[0])) + " bp in length.")
Next, we want to digest this scaffold into guides. This uses the al_diggesttarget function contained in eating.py to produce a generator of scores.
In [7]:
cutslist = al_digesttarget(longscaffold)
In this version of the script, the output from al_digesttarget isn't especially readable, but for reference:
In [8]:
[item for item in al_digesttarget([longscaffold[0][0:1500]])]
Out[8]:
Next, we'd like to take the 20mers extracted and score them against the entire Xenopus laevis genome. These lines score each guide variable region for specificity using the xl71 BLAST database and the xl71genomedict dict.
In [ ]:
def multiscore_pool(x):
score = al_scoreguide(x, "xl71", xl71genomedict)
return (score[0], score[1])
In [11]:
http://sebastianraschka.com/Articles/2014_multiprocessing_intro.html#An-introduction-to-parallel-programming-using-Python's-multiprocessing-module
pool = mp.Pool(processes=2)
results = [pool.apply(multiscore_pool, args=(x,)) for x in cutslist]
pickle.dump(results, open( "finalpicklescores.pkl", "wb" ))
pool.close()
The format of the resulting data is (score, guide).
In [53]:
results[0:20]
Out[53]:
The scores in this object are an ordered list, with all HpaII scores first, then all BfaI scores and finally all ScrFI scores. We are interested in the distribution of scores along the DNA fragment, irrespective of the enzyme used to generate them. Thus, we want to rearrange the list with all scores from 5' to 3'.
In [93]:
import copy
a = []
for (score, details) in results:
a.append(int(details.name)) # The guide's name attribute contains its position in bp
resultssorted = zip(results, a)
resultssorted = sorted(resultssorted, key=itemgetter(1), reverse=False)
resultssorted = [item for item, null in resultssorted]
In [125]:
resultssorted[:5]
Out[125]:
In [126]:
resultssorted[-5:]
Out[126]:
Let's extract the scores and plot their distribution on a histogram.
In [127]:
scores = [score for score, details in resultssorted]
In [128]:
def plot_score_histogram(scores):
'''
Input is a list of scores only (as ints)
'''
path = '/Library/Fonts/Microsoft/Arial.ttf'
prop = matplotlib.font_manager.FontProperties(fname=path)
matplotlib.rcParams['font.family'] = prop.get_name()
bins = range(0,106,5)
figure()
hist(scores, bins, color="gray")
tick_params(axis=u'both', labelsize=18)
#savefig('Scaffold score distribution.pdf', format="pdf")
In [130]:
plot_score_histogram(scores)
So, there are ~5000 guides that are quite non-specific (score <= 4) and >14,000 guides that have a score of 100 and a further 4000 that score between 95 and 99.
To make a library of useful guides, we'd like to PCR through continuous clusters of the highest-scoring ones. Our oligonucleotide vendor (IDT) has a minimum order of 288 oligos (3x 96-well plates) on a small and relatively inexpensive scale (5 pmol). To work within this limitation, we'd like to pick out 144 possible regions to PCR-amplify.
If we are only willing to accept guides with a score of 100, we'd predict that our 144 PCR products will be short (there are probably few long spans of perfect-scoring guides). However, if we relax our requirement to >=99, we may get longer PCR products and thus more guides in our library. How does this scale at different cutoffs/thresholds?
In [215]:
def find_clusters_by_cutoff(resultssorted, x):
starts=[]
ends=[]
previtemgood = 0
for index, (score, details) in enumerate(resultssorted):
if score >= x and previtemgood ==0 and len(details) >= 20: #this avoids guides that are shorter than 20 bp (from where an enzyme cuts twice in close proximity)
starts.append((index, score, int(details.name)))
previtemgood = 1
elif score >= x and previtemgood == 1 and len(details) >=20:
None
elif previtemgood == 1:
previtemgood =0
ends.append((index-1, resultssorted[index-1][0], int(resultssorted[index-1][1].name)))
run_positions = zip(starts, ends)
goodruns_length = sorted([end - start for (start, i, j), (end,l,m) in run_positions], reverse=True)
return (goodruns_length, run_positions)
In [216]:
threshold = range(0, 105, 5)
probeyield = []
for item in threshold:
probeyield.append((item, sum(find_clusters_by_cutoff(resultssorted, item)[0][0:143])))
In [217]:
print(probeyield)
In [148]:
%pylab inline
figure()
plot([b for b, c in probeyield], [c for b, c in probeyield], "o")
Out[148]:
Our "yield" of guides descends steadily from a cutoff of >=5 to a cutoff of >=95, then drops from 2894 guides produced at a cutoff of 95 to 1719 at 100. So, a cutoff of >=95 balances specificity and yield.
In [205]:
threshold = 95
runs = find_clusters_by_cutoff(resultssorted, threshold)[1]
#(countofguides, (startguidenumberfrom5', startscore, startpositionbp), (endguidenumberfrom5', endscore, endpositionbp))
goodruns = sorted([((i_end - i), (i, s, pos), (i_end, s_end, pos_end)) for (i, s, pos), (i_end, s_end, pos_end) in runs], reverse=True)
We next asked what happens if we concentrate the guides into a smaller region. To test this, we cut the input DNA into sections of 1/7 the ~21MB starting length and asked how many guides would be obtained if 144 PCR products were designed within each of those subregions.
In [ ]:
probeyield = []
x = 95
fraction = 7.0
overlap = 2.0
region_to_extract = len(resultssorted)/fraction
for i in [float(item)/overlap for item in range(int(overlap*fraction+2.0))]:
goodruns = find_clusters_by_cutoff(resultssorted[int(region_to_extract*i):int(region_to_extract*(i+1))], x)[0]
probeyield.append((i, int(region_to_extract*i), sum(goodruns[0:143])))
if sum(goodruns[0:143]) == 0:
break
In [248]:
probeyield
Out[248]:
The final 1/7 of the scaffold has the densest guide yield.
In [ ]:
#Modify resultssorted to only include the 3.4MB region used. (18121076 to (21505465+786) = 21506251)
resultssorted = [item for item in resultssorted if int(item[1].name) >= 18121076 and int(item[1].name) <= 21506251]
scores = [score for score, details in resultssorted]
We want to make sure we design primers to make PCR products spanning as many guides as possible. The challenge is making sure we prime in a way that covers all of the high-specificity guides and none of the adjacent low specificity guides. If a "good" and a "bad" guide are very close together, the constraints on where primers can be designed against are tight.
In [ ]:
# Set up the input for primer3:
# Sequence available to PCR:
guide_count = []
amps_in_3MB = []
for index, item in enumerate(goodruns[0:400]):
left_outside = item[0][0].id[-1]
left_inside = item[2][0][1].id[-1]
if left_outside == "F" and left_inside == "R":
permissible_start = int(item[0][0].name) + 10
required_start_absolute = int(item[2][0][1].name) +14
elif left_outside == "R" and left_inside == "R":
permissible_start = int(item[0][0].name) + 1
required_start_absolute = int(item[2][0][1].name) +14
elif left_outside == "R" and left_inside == "F":
permissible_start = int(item[0][0].name) + 1
required_start_absolute = int(item[2][0][1].name) +18
elif left_outside == "F" and left_inside == "F":
permissible_start = int(item[0][0].name) + 10
required_start_absolute = int(item[2][0][1].name) +18
else:
print("error on left")
right_inside = item[2][-1][1].id[-1]
right_outside = item[0][1].id[-1]
if right_outside == "F" and right_inside == "R":
permissible_end = int(item[0][1].name) + 19
required_end_absolute = int(item[2][-1][1].name) + 2
elif right_outside == "R" and right_inside == "F":
permissible_end = int(item[0][1].name) + 10
required_end_absolute = int(item[2][-1][1].name) + 8
elif right_outside == "R" and right_inside == "R":
permissible_end = int(item[0][1].name) + 10
required_end_absolute = int(item[2][-1][1].name) + 2
elif right_outside == "F" and right_inside == "F":
permissible_end = int(item[0][1].name) + 19
required_end_absolute = int(item[2][-1][1].name) + 8
else:
print("error on right")
amp = longscaffold[0][permissible_start:permissible_end]
# Bounds that need to be included in PCR product :
required_start_relative = required_start_absolute-permissible_start
required_end_relative = required_end_absolute - permissible_start
amp.dbxrefs=((required_start_relative, required_end_relative))
# Set up some other stuff:
amp.name =str(item[0][0].name)
amp.id =str(item[0][0].name)
amp.description=str(item[1])
amp.seq.alphabet = IUPACAmbiguousDNA()
if "NNNNN" in amp.seq: # Exclude if it has runs of Ns
None
#print amp.name + " contains ns " + str(item[1])
else:
amps_in_3MB.append(amp)
guide_count.append(item[1])
amps_in_3MB_gen = (i for i in amps_in_3MB)
print sum(guide_count[0:144])
In [ ]:
with open("primerlist.txt", "w") as primerlist:
primerlist.write("Sequence_id\tforward_seq\tforward_start\tforward_length\tforward_tm\tforward_gc\treverse_seq\treverse_start\treverse_length\treverse_tm\treverse_gc\tinput_seq_length\tPCR_product_length\tGuides_Contained\n")
primerlist.close()
for item in amps_in_3MB:
current_amp = item
primerdict = al_primersearch(current_amp)
al_collect_good_primers(item, primerdict)
This outputs to "primerlist.txt" the final list of primers to PCR from the Xenopus genome, listed in Table 3 (this is tab-separated so can be imported into Excel). After PCR, the products are pooled and subjected to the EATING protocol.